Adaptive Mesh Refinement Computation of Solidification Microstructures using 

Dynamic Data Structures 



Nikolas Provatas^'^, Nigel Goldenfeld^ , and Jonathan Dantzig^ 
^University of Illinois at Urhana- Champaign, Department of Physics 
1110 West Green Street, Urbana, IL, 61801 
^ University of Illinois at Urbana-Champaign, Department of Mechanical and Industrial Engineering 

1206 West Green Street, Urbana, IL, 61801 
(February 1, 2008) 

and speed. Experiments on dendrite evolution by Glicks- 

man and coworkers on succinonitrile (SCN) and 
more recently Pivalic Acid (PVM) as well as other 
transparent analogues of metals have provided tests of 
theories of dendritic growth, and have stimulated consid- 
erable theoretical progress . These experiments have 
clearly demonstrated that in certain parameter ranges 
the physics of the dendrite tip can be characterized by a 
steady value for the dendrite tip velocity, radius of curva- 
ture and shape. Away from the tip the time-dependent 
dendrite exhibits the characteristic sidebranching as it 
propagates. 

The earliest theories of dendritic growth solved for the 
diffusion field around a self-similar body of revolution 
propagating at constant speed ■ In these studies the 
diffusion field is found to determine the product of the 
dendrite velocity and tip radius, but neither quantity by 
itself. Adding capillarity effects to the theory predicts a 
unique maximum growth speed j9j. Experiments, how- 
ever, have shown that these theories do not represent the 
correct operating state for real dendrites. 

The introduction of local models of solidification 
brought further insight to the steady state dendrite prob- 
lem [pO|-p^. These models describe the evolution of the 
interface, incorporating the physics of the bulk phases 
into the governing equation of motion of the interface. 
A remarkable breakthrough of these models was to show 
that a steady-state dendrite velocity is obtained only if 
a source of anisotropy (e.g., in the interfacial energy) is 
present during dendritic evolution. The dendrite steady- 
state tip velocities appear in a discrete rather than con- 
tinuous spectrum of values, making the role of anisotropy 
of great importance in the description of the dendrite 
problem, both in the local models and the full moving 
boundary problem P, p"4yT^ . It was further shown that 
only the fastest of a spectrum of steady state velocities 
is stable, thus forming the operating state of the den- 
drite. This body of theoretical work is generally known 
as microscopic solvability theory. 

Dendritic sidebranching is widely believed to be caused 
by thermal fluctuations, which enter solidification models 
in the form of random noise possessing specific features 
p6|-p3. However, the precise mechanism of sidebranch 



We study the evolution of solidification microstructures 
using a phase-field model computed on an adaptive, finite 
element grid. We discuss the details of our algorithm and 
show that it greatly reduces the computational cost of solv- 
ing the phase-field model at low undercooling. In particular 
we show that the computational complexity of solving any 
phase-boundary problem scales with the interface arclength 
when using an adapting mesh. Moreover, the use of dynamic 
data structures allows us to simulate system sizes correspond- 
ing to experimental conditions, which would otherwise require 
lattices greater that 2^^ x 2^^ elements. We examine the con- 
vergence properties of our algorithm. We also present two 
dimensional, time-dependent calculations of dendritic evolu- 
tion, with and without surface tension anisotropy. We bench- 
mark our results for dendritic growth with microscopic solv- 
ability theory, finding them to be in good agreement with 
theory for high undercoolings. At low undercooling, however, 
we obtain higher values of velocity than solvability theory at 
low undercooling, where transients dominate, in accord with 
a heuristic criterion which we derive. 



I. INTRODUCTION 

Modeling solidification microstructures has become an 
area of intense study in recent years. The properties of 
large scale cast products, ranging from automobile engine 
blocks to aircraft components and other industrial appli- 
cations, are strongly dependent on the physics that occur 
at the mesoscopic and microscopic length scales during 
solidification. The main ingredient of the solidification 
microstructure is the dendrite, a snowflake-like pattern 
of solid around which solidification proceeds. The micro- 
scopic properties of such cast products are determined 
by the length scales of these dendrites, and for this rea- 
son understanding the mechanisms for pattern selection 
in dendritic growth has attracted a great deal of inter- 
est from the experimental and theoretical community. In 
particular, a great deal of research has been undertaken 
to understand such issues as dendrite morphology, shape 
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amplification does not seem to be fully understood. It 
would appear that the manner in which thermal noise is 
amplified may depend on the overall dendrite morphol- 
ogy. For instance, it was shown that noisy fluctuations 
traveling down a paraboloid of revolution do not pro- 
duce sidebrach amplitudes consistent with experiments 
I p^ , while fluctuations traveling down an initially missile- 
shaped dendrite amplify into sidebranches comparable to 
some experiments Karma also investigated the ad- 
dition of interface fluctuations However, this source 
of noise only becomes relevant at high velocities. 

The theoretical foundation around which most theo- 
ries of solidiflcation are based is the time-dependent Ste- 
fan problem. This theory describes the evolution of the 
thermal or solutal diffusion fleld around the solidifica- 
tion front, along with two accompanying boundary condi- 
tions. The first boundary condition relates the velocity of 
the moving front to the difference in thermal fiuxes across 
the solid-liquid interface. The second, called the Gibbs- 
Thomson condition, relates the interfacial temperature 
to the the thermodynamic equilibrium temperature, the 
local interfacial curvature and interface kinetics. The in- 
terface kinetics term adds a non-equilibrium correction 
to the interface temperature, usually assumed to be in 
local equilibrium for a given curvature. Solving the Ste- 
fan problem numerically has traditionally involved front 
tracking and lattice deformation to contain the interface 
at predefined locations on the grid 20|. This method 
is generally complicated to implement accurately and re- 
quires much effort. Moreover, it can be inefficient in han- 
dling coalescence of two or more interfaces. 

The solution of the Stefan problem has been made 
more tractable with the introduction of the phase-field 
model. The phase-field model avoids this problem of 
front tracking by introducing an auxiliary continuous or- 
der parameter (j){r) that couples to the evolution of the 
thermal or solutal field. The phase field interpolates be- 
tween the solid and liquid phases, attaining two different 
constant values in either phase, (e.g., ±1) with a rapid 
transition region in the vicinity of the solidification front. 
The level set of (?!)(r) = is identified with the solidifica- 
tion front, and the subsequent dynamics of cj) are designed 
to follow the evolving solidification front in a manner that 
reproduces the Stefan problem pl]-^. 

The price to be paid for the convenience of the or- 
der parameter is the introduction of a new length scale 
W which represents a boundary layer over which the or- 
der parameter changes sign. This distance is referred to 
as the interface width, and does not appear in the Ste- 
fan problem. As such, one requirement of the phase-field 
model is to recover the Stefan limit in a manner that is in- 
dependent of the interface width as W approaches some 
appropriate limit, considerable work has been done to 
relate W to various parameters of the phase-field model 
in order to establish a mappin g b etween the phase-field 
model and the Stefan problem ||2^,^,^ . While the for- 
mal nature of these mappings does not seem to be very 
sensitive to the precise form of the phase- field model [p3| , 



different asymptotic limits of the phase-field parameters 
can lead to widely varying complexity in the numerical 
implementation. 

The introduction of the interface width W makes the 
phase-field model prohibitively expensive to simulate nu- 
merically for large systems, since the grid spacing must 
be small enough everywhere that the phase-field model 
converges to the the sharp interface limit |2^,^. Cagi- 
nalp and Chen |Q showed rigorously that the phase 
field model converges to the sharp interface limit when 
the interface width (and hence the grid spacing) is much 
smaller than the capillary length. While this result is nec- 
essary to establish that the phase-field model does map 
onto the Stefan problem, the parameter values required 
to realize the asymptotic limit can be computationally 
intractable. Experimentally, the physical sizes required 
to contain realistic microstructures can be many times 
the size of the thermal diffusion length, which in turn 
can be orders of magnitude greater than W Thus, since 
Axinin < W, computing in the hmit of a ^ does 
not allow one to simulate experimental systems. 

Recently Karma and Rappel presented a different 
asymptotic analysis performed in powers of the ratio of 
the interface width to the diffusion length a/Vn, taken to 
be equal in both phases. Their procedure offers two com- 
putational advantages. The flrst is that is allows one to 
simulate the phase-fleld model with zero interface kinet- 
ics, without the need to make W ^ 0. Specifically, this 
limit, as well as a non-zero kinetics limit, can be simu- 
lated with an interface width W (and hence the grid spac- 
ing) of order the capillary length, a much more tractable 
regime. Simulating solidification microstructures in the 
limit of zero interface kinetics is important because most 
experiments performed at low undercooling in materials 
such as succinonitrile are in this limit Karma and 
Rappel tested their asymptotics by comparing their sim- 
ulations to the results of microscopic solvability theory, 
finding excellent agreement down to dimensionless un- 
dercoolings as low as 0.25. 

A recent extension of Karma and Rappel's analysis by 
Almgren [^5| also promises to allow similar asymptotics 
to be performed on a two-sided model of solidification 
p5| , i.e. when the diffusivities in the solid and liquid 
differ, relevant in the study of directional solidification of 
binary mixtures. 

The theory of level sets has also recently re- 

emerged as another effective tool that shows great po- 
tential in modeling dendritic growth. While related to 
the phase-field model, level-set theory does not require 
the presence of a thin interfacial with W, thus greatly 
reducing the stringent grid requirements posed by con- 
ventional phase-field models. To date, however, level set 
methods have not been benchmarked with solvability the- 
ory or other theoretical prediction for Stefan problems. 

While expanding the horizon of solidification model- 
ing, phase-field modeling has still been limited to small 
systems sizes, even when solved by adaptive algorithms 
|p8[. The main problem is the presence of an interface 
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region with a minimal length scale that must be resolved. 
For a typical microstructures grown at dimensionless un- 
dercooling of 0.1 or less, the ratio of the system size to 
this minimal grid spacing can be greater that 2^^. With 
this restriction most numerical methods will naturally 
fail. What is needed to go beyond this limitation is an 
effective adaptive technique |38-4^ which dynamically 
coarsens the grid spacing away from the front. 

In this paper we present a new, computationally effi- 
cient adaptive-grid algorithm for solving a class of phase- 
field models suitable for the study of phase-boundary 
evolution. We study dendritic solidification modeled us- 
ing two coupled fields, one for the order parameter and 
the other for the thermal field. Our algorithm effectively 
combines and implements ideas of adaptive-mesh refine- 
ment with ideas of dynamic data structures, allowing us 
to enlarge the window of large-scale solidification model- 
ing. 

The outline of this paper is as follows: In section one we 
introduce the physical model to be examined, summariz- 
ing its properties and its various limits. In section two we 
present a detailed description of our algorithm. In section 
three we present results on CPU-scalability of our algo- 
rithm and examine issues of grid convergence and grid 
anisotropy on our solutions. In section four we present 
results of dendritic growth with and without the presence 
of anisotropy in the surface energy. We show that for high 
undercooling, dendrites grown with our method converge 
to tip speeds in agreement with microscopic solvability 
theory. At low undercooling, however, we do not find 
agreement with steady state solvability theory, owing to 
long-lived transients in the thermal field evolution. In 
section four we conclude and discuss our results. 



II. THE PHASE-FIELD MODEL 

We model solidification using a standard form of phase- 
field equations which couple a thermal field to an or- 
der parameter field (p via a double well potential |^,^ . 
We begin by rescaling the temperature field T hy U = 
cp{T — Tm)/L, where cp is the specific heat at constant 
pressure, L is the latent heat of fusion and Tm is the 
melting temperature. The order parameter is defined by 
0, where we define (f> — 1 in the solid phase and (f> — —1 
in the liquid phase. The interface is defined by = 0. 
We rescale time throughout by Tq, a time characteriz- 
ing atomic movement in the interface. Length is rescaled 
by Wo, a length characterizing the liquid-solid interface. 
With these definitions, the model is written as 



Itt 



2'dt 



(1) 



AHn)^^V- {A' (n) V0) + g (0) - XUP' (0) 



where D = aTo/W^ and a is the thermal diffusivity. The 
function /(</>, U; A) = 5'(</>) — XUP'{<j)) is the derivative of 
the double-well potential with respect to (f> and couples 
the U and (j) fields via the constant A. The primes on the 
functions g{(j)) and P{(j)) denote derivatives with respect 
to (f). Following Karma and Rappel |3^, anisotropy has 
been introduced in Eqs. (|l|) by defining the width of the 
interface to be W{n) = WoA{n) and the characteristic 
time by T{n) — ToA'^{n), with A{n) £ [0, 1] given by 



A{n) = (1 - 3e) 
The vector 



4e 



l-3e |V0|^ 



(2) 



(3) 



defines the normal to the contours of the field, where 
and (^^y are defined as the partial derivatives of </> 
with respect to x and y. The variable e parameterizes the 
deviation of W{n) from Wo and represents the anisotropy 
in the interface energy of the system. We note that this 
definition of anisotropy is not unique |Q , but we expect 
results to be similar for other definitions of anisotropy. 

In simulating the phase-field model we adopt the point 
of view that the order parameter field is a compu- 
tational tool whose main purpose is to eliminate front 
tracking. As such we would like to simulate the model 
given by Eqs. (|l|) with Wo as large as possible. At the 
same time we would like the behavior of the model out- 
side the boundary layer defined by to describe the Ste- 
fan problem as closely as possible. To this end, we relate 
the parameters of the phase-field model according to Ref. 
p2[ , valid in the asymptotic limit Wo <C a/Vc, where 
a/Vc is the diffusion length and Vc is a characteristic ve- 
locity of the front defined by 0. 

The specific asymptotic limit we model is one where 
the [/-field satisfies 



dU_ 
'dt 



(4) 



everywhere away from the interface, while at the inter- 
face, the gradient of U satisfies 




dU 




lilt 





(5) 



where Vint is the velocity normal to the interface, denoted 
by Xint- The notation ± denotes the solid/hquid side of 
the interface, respectively. The description of the Stefan 
problem is completed by the Gibbs-Thomson condition 
and the interface kinetics condition 



C^(^int) = -d{n)K - f3{n)Vn, 



(6) 



where d{n) is the capillary length, k is the local curvature 
and (3{n) is the interface attachment kinetic coefficient. 
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all assumed to be in dimensionless form according to the 
above definitions. The capillary length is related to the 
parameters of Eqs. (m) by 



d{n) 



ai- 



Wo 

X 



[A{n) + d'oA{n)] 



(7) 



where ai = 0.8839 for the particular form of the free 
energy defined in Eqs. (0) and 9 is the angle between 
n and the x-axis. The kinetic coefficient is given by 



XWo 



1 



Aa2 



(8) 



where 02 = 0.6267 for our choice of the free energy func- 
tional citeKar95. One remarkable feature of Eqs. (j^) and 
is that Wo, To and A an be chosen to simulate arbi- 
trary values of /3, for Wo of order do ■ In particular, setting 
A = D/a2 allows us to compute the phase-field model in 
the limit of the Stefan problem where /3 = 0. This 
is also an appropriate value for SCN, especially at low 
undercooling. 

Equations (0) and for /? and do can be related to 
a wide class of free energies via the parameters ai and 
a2 |3^, which are related to integrals that depend on 
g{4>o), P{,4'o) and d(j)o/dx, where 0o is the lowest order 
description of the order parameter field (fi and and is a 
solution of the equation 



d'^(l)o _ dg{(i)o) 

dx^ d(t>n 



0. 



(9) 



We also note that these asymptotics are a special case 
of a more general asymptotic analysis performed by Alm- 
gren , which relates the parameters of the phase- field 
model to those of the Stefan problem in the case of un- 
equal diffusivities in the solid and liquid phases. In this 
case, the asymptotics provides an additional set of con- 
straints on the admissible functions P'{4>), g'{<p)i aiid 
hence ai and 02- 



III. THE ADAPTIVE-GRID ALGORITHM 

The main computational challenge of simulating 
Eqs. (|^) involves resolving two competing length scales: 
the lattice spacing dx on which the simulation in per- 
formed and the physical size of the system L b ■ Even with 
improved asymptotics, dx must remain relatively small, 
while Lb must be extremely large in order to make pos- 
sible computations of extended solidification microstruc- 
tures. Moreover, the main physics of solidification (and 
the evolution of most phase-boundary problems) occurs 
around an interface whose area is much smaller than the 
full computational domain. Near this interface the or- 
der parameter varies significantly, while away from the 
interface variations in are small. Meanwhile, the ther- 
mal field U extends well beyond the interface and has 
much more gradual variation in its gradients, permitting 



a much coarser grid to be used to resolve U. The most 
obvious manner to overcome this problem is to use a 
method that places a high density of grid points where 
the interface of </) or [/ varies most rapidly and a much 
lower grid density in other regions. Furthermore, the 
method must dynamically adapt the grid to follow the 
evolving interface |38|-|4l|, while at the same time main- 
taining a certain level of solution quality. 

We solve Eqs. (|l|) using the Galerkin finite element 
method on dynamically adapting grids of linear, isopara- 
metric quadrilateral and triangular elements. The grid 
is adapted dynamically based on an error estimator that 
utilizes information from both the and U fields. We 
wrote our code in FORTRAN 90 (F90), taking advan- 
tage of the efficiency of FORTRAN 77 while using ad- 
vanced C-like features, such as data structures, derived 
data types, pointers, dynamic memory allocation and 
modular design to conveniently adapt the grid and the 
solution fields. 

In the broadest sense, our algorithm performs func- 
tions that can be divided into two classes. The first deals 
with the establishment, maintenance and updating of the 
finite element grids, and the second with evolving and 
U on these grid, according to Eqs. (|l|). We presently de- 
scribe these classes, the adaptive grids, the finite element 
procedure, and the interconnections of these processes. 



A. The Finite Element Grids 

The first class of functions in our algorithm centers 
around maintaining a grid of finite elements on a data 
structure known as a quadtree p^-^^ which replaces the 
standard concept of a uniform grid as a way of repre- 
senting the simulational grid. The quadtree is a tree- 
like structure with branches up to a prespecified level. 
Branches of the quadtree are themselves data structures 
that contain information analogous to their parent, from 
which they branched, but one level down. Fig. |l| illus- 
trates the structure of a quadtree as well as the relation 
between elements at different levels of refinement. Every 
entry on the quadtree contains information pertaining 
to a 4-noded isoparametric quadrilateral finite element. 
This information includes the following: 

• values of (p and U at the four nodes 

• the nodal coordinates of the element 

• the level of refinement of the element on the quadtree 

• the value of the current error estimate 

• The element number, which contains information about 
the coordinates of the element and its level of refinement 

• an array mapping the element's four nodes onto the 
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entries of a global solution array 



ing information: 



• pointers to the element's nearest neighbors sharing a 
common edge and at the same level of grid refinement 

• a variable that determines whether or not an element 
contains further sub-elements which we term child ele- 
ments 

• pointers to an element's child elements 

• a pointer to the parent element from which an element 
originates 

The elements of the quadtree can be refined by split- 
ting into four child elements, each sharing the same par- 
ent element one level down on the quadtree and each 
with its own unique information, as outlined above. A 
parent element and it's four child elements are referred 
to as a family. Refinement produces a finer mesh within 
the confines of the original parent grid by bisecting each 
side, as shown in Fig. (ij). Unrefinement, which consists 
of fusing the four child elements back into the parent, 
has the opposite effect, locally creating a coarser mesh. 
Both refinement and unrefinement proceed via dynamic 
memory allocation, making our code scalable. We note 
that unrefinement can occur only if the child elements 
do not possess further child of their own. Also, in or- 
der to avoid having regions of very different refinement 
bordering each other, we impose the restriction that any 
two neighboring quadrilateral elements may be separated 
by no more than one level of refinement (see Fig. (0)). 
We define the level of refinement of an element by le such 
that a uniform grid at a refinement level le would contain 
2'= X 2'= grid points in a physical domain Lb x Lb- 

Special cases where an element has no children, a miss- 
ing neighbor, or no parent are handled by null pointers. 
The latter case occurs only for the root of the quadtree. 

All elements at a given level of refinement on the 
quadtree are "strung" together by a linked-list of point- 
ers, referred to as the G-lists. There are as many G-lists 
as there are levels of refinement in the quadtree. Each 
pointer in the G-list points to (accesses) the location in 
memory assigned to one element of the quadtree. The 
purpose of the G-list is to allow traversal of the quadri- 
lateral elements sequentially by level, rather than by re- 
cursively traversing quadtree from the root down, a pro- 
cedure which is memory intensive and relatively slow. 

Alongside the main grid of quadtree elements, the code 
maintains two independent grids representing special lin- 
ear isoparametric triangular and rectangular elements. 
These elements are used to connect the extra nodes that 
arise when two or more quadrilateral elements of differing 
refinement levels border each other. These element types 
are referred to as bridging elements. They are maintained 
as two linked- lists of derived data types, one containing 
information about triangular elements and the other rect- 
angular. Elements of both these grids include the foUow- 



• the values of (p and T at the three nodes (four for rect- 
angles) of the element 

• the nodal coordinates 

• node numbers that map the element's nodes onto the 
global solution array 

The types of bridging triangles and rectangles that can 
occur are enumerable and shown in Fig. (H). 

The main set of operations performed on the grids de- 
scribed above concern refinement of the finite element 
mesh as a whole. The refinement process is performed 
only on the quadrilateral mesh. The triangular and rect- 
angular grids are established after this process is com- 
pleted (see Fig. (|l|)). To refine the grid the code traverses 
the elements of the quadtree, refining (unrefining) any 
element whose error function, discussed below, is above 
(below) a critical value ah{(Ji). We also note that fusion 
of four quadrilateral elements can occur only if all four of 
its children's error functions are below the critical value 
ai, where ai < ah- We found that if ai = ah the grid sets 
into oscillations, where large numbers of elements become 
alternatively refined at one time step, then unrefined at 
the next. 

The processes described thus far are grouped into mod- 
ules that encapsulate various related tasks, and which 
can cross-reference each other's data and instructions. 
The module highest up in the hierarchy contains the def- 
inition of the quadtree data structure and routines that 
construct the initial uniform grid, refine and unrefine in- 
dividual quadrilateral elements, and set the initial condi- 
tions. Another module constructs the G-lists. It contains 
routines that construct the initial G-list from initial uni- 
form quadtree data structure, as well as add or delete 
element pointers from the G-list as elements are created 
or deleted from the quadtree. Another module accessing 
both the previous ones' data structures has the role of 
creating the triangular and rectangular element grids. It 
contains definitions for creating triangular and rectangu- 
lar elements data structures and routines that search the 
quadtree, building the linked lists of triangles and rect- 
angles that make up these grids. The main program is 
contained in its own module and contains the driver pro- 
gram that creates the initial grids, G-lists and triangular 
and rectangular element types. The driver program also 
sets into motion the final link in the simulation, which 
evolves (/) and U and periodically adapts the dynamic 
grid by calling procedures described above. A flowchart 
of these processes is shown in Fig. (S). 
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B. The Finite Element Formulation 

The integration of Eqs. is done by the final module 
in the code. This module performs four main processes: 



1. Maps the internal element node numbers to the 
indices of a global solution vector. The ip-&e\d is 
mapped onto the odd numbers, (1, 3, 5, . . .), while 
U is stored on the even numbers of the global solu- 
tion vector (2, 4,6,.. .) 

2. Advances the U and cj) field- vectors by Nr time 
steps on the finite element grids defined above 

3. Calculates an error function for each element of the 
quadtree, based on error estimate of the quadrilat- 
eral elements 



for i = 1,2,5, N , where represents the area of 
an element e. Substituting Eqs. into Eqs. (p^), we 
obtain two linear algebraic equations for (pi and Ui, i ^ 
1, 2, 3, . . . , iV in the element e. 

We next define {^Y — </'2, ^3, • ' ' j ^w)^ and 
{U}'^ — {Ui, U2, U3, ■ ■ ■ , Un)'^ , where the superscript T 
denotes transpose, making {^}'^ and {U}^ column vec- 
tors. The linear algebraic statement of the finite element 
form of Eqs. (|l]) then becomes 

[cmv)^ = im + m + {f-, xy (13) 

where the matrices [C], [C], [A], [M] and [E] and the 
vector {F; A}^ are given by 



[C] = [NY [N]dxdy, 



(14) 



4. Invokes routines in the modules described above to 
refine the grid according to this error estimator 



[C] = / [NY [N]A\e{cj,^))dxdy, 

J fie 



(15) 



Steps (l)-(4) are repeated until a sufficient time evolu- 
tion of the microstructure is established. The variable N^ 
is set such that the interface remains within the regions of 
fine mesh between regriddings, which we typically choose 
to be 100 time steps. Step (1) involves searching all el- 
ements, and their neighbors, and assigning each node a 
unique number, that will have a counterpart on a global 
solution vector. 

The finite element discretization of Eqs. ([l|) is done 
using Galerkin's weighted residual method p5[. The 
method begins by assuming that <j) and U are interpo- 
lated within an element as 



N 

E 

i=l 



N 

E 

i=l 



U^N,ix,y) (10) 



where (fif and Uf are the field values at the N nodes of 
the element e, and their interpolated values in its interior. 
The functions Ni{x,y) are standard linear interpolation 
functions appropriate to the element being used [ 46| , and 
satisfy 



Nt{xj,yj) = 5, 



(11) 



where Si^ is the Kroneker delta. Rewriting the differen- 



tial equations for cj) in Eqs. (|1|) as L^cf) = 0, and of the 
?7-equation as LjjU — 0, the Galerkin method requires 
that 



Ni{x,y)L^(t)''{x,y)dxdy = 



/ N,{x, y)LuU''{x, y)dxdy = 0, 
Jo, 



(12) 



[^] = - / {[NV[N.] + [^^[Ny]) dxdy, (16) 



[M] = - / {[Nf[N,] + [Nf[Ny]) A\e{r))dxdy, 



(17) 



[E] = - {[Nf[N,] - [Nf[Ny]) A{emM0m)dxdy, 



(18) 



{F;Xr^f [Nff{cP',U-;X)dxdy, (19) 

where [N^] , [Ny] denote the partial derivatives of the vec- 
tor of shape functions with respect to x and y, respec- 
tively. The function A is just Eq. (|^) rewritten in terms 
of the angle 9 that the normal to the interface makes with 
the X-axis. Specifically, defining 



tan 6* 

then 

A{9m) = (l-3e) 



4e (l-f tan^ 



1 - 3e (1 tan 



(20) 



(21) 



while Lo{9) is proportional to the derivative of A{9), and 
is given by 
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w(6'(0")) = 16e 



tan 0(1 — tan^ 
(1 + tan2 0): 



(22) 



We use a lumped formulation for the matrices [C] and 
[C] In this procedure, the row vector of shape func- 
tions, [N] in Eq. (14) is replaced by the identity row 
vector [I] — [1,1,1,---]. The resuhing matrix [C] then 
consists of identical columns, each of which contains the 
element Ni{x, y) in the position of the i"^ row. A lumped 
term is then defined as a diagonal matrix whose entries 
take on the value 



^ nodes „ 

= T / N,{x,y)dxdy. 



4 ^ 

i=l 



(23) 



The use of a lumped matrix for [C] allows us to assemble 
a diagonal matrix for the left hand side Eqs. (|l3|), stored 
as a one-dimensional vector rather than two-dimensional 
matrices that would be required if we used the consistent 
formulation for the assembly of the [C] matrices. In- 
deed, microstructures evolving at low undercooling can 
produce interfaces with over 2 x 10^ elements, making 
the storing of 2 x 10^" matrices impossible. 

The global {</)} (obtained after assembly of the element 
equations in field in Eqs. (13)) is time-stepped using us- 
ing a forward difference (explicit) time scheme. For each 
time step of the (j) field, the global U field is then solved it- 
eratively using a Crank-Nicholson scheme. Convergence 
of {U}n+i is obtained in a few iterations. 



C. The Error Estimator 

Regridding is based on an error estimator func- 
tion, which is obtained following Zienkiewicz and Zhu 
p6| , based on the differences between calculated and 
smoothed gradients of the and U fields. Specifically, 
we define the composite field 



* = + 7C/ 



(24) 



where 7 is a constant. We discuss the selection of 7 in 
more detail below. This definition allows us to regrid ac- 
cording the requirements of both the (f> and U field, as 
opposed to using only gradients of the 0-field in estab- 
lishing the grid ||3^. Since it is and U that are being 
calculated, and not their gradients, we do not expect the 
gradient of ^' to be continuous across element boundaries, 
due to the order of the interpolation used. Thus we ex- 
pect the difference between the calculated and smoothed 
(continuous across element boundaries) gradients to pro- 
vide a reasonable estimate of error. This method ap- 
propriately meshes regions of both steep gradients and 
regions where the and U fields change rapidly. 
We define the error estimator function e as 



(25) 



where qc and Qs are the calculated and smoothed gra- 
dients of ^E" respectively. The smoothed gradients are 
calculated to be continuous across element boundaries. 
To determine qg we assume it to be interpolated in the 
same way as the (p and U fields, namely 



qs = [N]{Q'} 



(26) 



where [N] is the row vector of element shape functions, 
and {Q*} is a 4 X 2 matrix whose columns represent the 
nodal values of fluxes of 5* in the x and y direction, re- 
spectively. To find {Q"} we use Galerkin's method, min- 
imizing the weighted residual 



[NYedn,= [NY {[N]{Q''} ~ q^)dn = (27) 



The calculation is simplified by lumping the left hand 
side of Eq. (|2^) , leading to 



[Nf[i]dn {Q^ 



[Nfqcdn, 



(28) 



where [1] = [1, 1, 1, • • • , TV]. Assembling Eq. (||) for 
all quadrilateral elements yields an equation for the 
smoothed gradients {Q}^ of the global field at all 
element nodes, of the form 



[D]{QV 



(29) 



where [D] is a diagonal matrix, due to "mass" lumping, 
and {Q}^ is a x 2 matrix for the global, smoothed flux. 

For the actual error updating on the elements of the 
quadtree we used the normalized error defined by 



(30) 



The domain of integration f2 in the denominator denotes 
the entire domain of the problem. Thus gives the con- 
tribution of the local element error relative to the total 
error calculated over the entire grid. 

Fig. (^ shows a snapshot at 10^ time steps into the 
a simulation of a thermal dendrite computed with our 
algorithm. The figure shows cj) and U as well as the cur- 
rent grid. The dendrite is four-fold symmetric, grown in a 
quarter-infinite space, initiated by a small quarter disk of 
radius Ro centered at the origin. The order parameter is 
defined on an initially uniform grid to be its equilibrium 
value (j)o{x) = — tanh((|2;| — Ro)/^/2) along the interface. 
The initial temperature decays exponentially from U = 
at the interface to — A as a? — *■ 00. The parameters set 
for this simulation are A = 0.70, D = 2, dt = 0.016 and 
A chosen to simulate /3 = 0. The system size is 800 x 800, 
with Axinin = 0.4, and about half of the computational 
domain in each direction is shown. Sidebranching is evi- 
dent, and arises due to numerical noise. This simulation 
was completed in approximately 15 cpu-hours on a Sun 
UltraSPARC 2200 workstation. 
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IV. SCALABILITY AND CONVERGENCE 
PROPERTIES OF THE ADAPTIVE-GRID 
ALGORITHM 

In this section we present results that iUustrate the 
convergence properties of solutions of Eqs. (P computed 
with our algorithm, the effect of grid-induced anisotropy 
of the adapting mesh, and the speed increase obtained 
by using an adapting grid. 



cpu time scales with arclength in the problem being con- 
sidered. We note that when interface convolutions be- 
come of order A ~ Axmim fine-grid regions separated 
by less than A will merge and the number of elements 
will stop growing locally. This makes the simulation of 
fractal-like patterns feasible as the arclength of the inter- 
face is bounded from above by x Lb- Finally, we note 
that adaptive gridding would especially improve the cpu 
performance of problems similar to spinodal decomposi- 
tion, where the total interface decreases with time. 



A. CPU-Performance 

We examined the cpu-scalability of our algorithm as a 
function of system size by growing dendrites in systems 
of various linear dimension Lb and measuring the cpu 
time required for the dendrite branches to traverse the 
entire system. Fig. || shows a plot of these data for a 
dendrite grown at undercooling A = 0.55 using the same 
parameters as in Fig. IJ. The minimum grid spacing has 
been set to Axmin = 0.4 in this data. Fig. ||, clearly 
shows that i?" L\. This relationship can be obtained 
analytically as follows. 

The number of calculations performed, per simulation 
time step, is proportional to the number of elements in 
the grid. This relationship is set in turn by the arclength 
of the interface being simulated multiplied by the diffu- 
sion length D/Vn- This product defines the arclength 
over which the highest level of refinement occurs. For a 
needle- like dendrite, the arclength is approximately Lb- 
Moreover, since the dendrite tip moves at a constant ve- 
locity Vn, then 



L 



B5 



(31) 



where i?" is a constant that depends on the details of the 
implementation of the algorithm used to evolve Eqs. (|^). 
The cpu time needed to compute the traversal time on a 
uniform grid, R", is found, by the same analysis, to be 



KAa 



r3 



(32) 



where i?" also depends on the implementation but is 
likely to be smaller than R"^. Thus, comparing our 
method with simulation on a uniform grid we obtain 



lim R'^jR'i = 
Lb^oo Lb 



(33) 



For larger systems, the adaptive scheme should always 
provide faster CPU performance regardless of implemen- 
tation. Indeed, any method that uses a uniform grid of 
any sort, will eventually be limited by memory require- 
ments as Lb becomes large. The arguments leading to 
Eq. (^Tj) can also be generalized to any problem of evolv- 
ing phase boundaries, always yielding the conclusion that 



B. Induced Lattice Anisotropy 

We tested the effective anisotropy of our dynamically 
adapting lattice in two independent ways. The first fol- 
lows the method outlined by Karma [Q. We fix the 
temperature far from the interface to be constant Too ev- 
erywhere, initially setting it to a critical value at which 
the isotropic surface energy just balances the bulk free en- 
ergy. For a specified background temperature, the crystal 
will only grow if its radius is greater than a critical value 
Ro- The radius Ro can be related to the background 
temperature through the total Gibbs free energy of the 
system, given by 



AG 



,LAT 



(34) 



where L is the latent heat of fusion, AT = T,n — Too, 
where Tm is the melting temperature. Too is the temper- 
ature far away from the interface, and a is the surface 
tension. Minimizing AG with respect to r yields Ro as a 
function of 6T as 



R* = do/ AT, 



(35) 



where do is the capillary length defined as do — 2(jTm/L- 
One finds an equilibrium shape of the interface when 
the background temperature field AT (written in terms 
of U) is adjusted dynamically so as to maintain the veloc- 
ity of the interface at zero as measured long the a;-axis. 
Thus, AT is increased if the velocity decreases, and de- 
creased if it grows. The effective anisotropy is inferred 
by fitting the computed interface to an equation of the 
form 



R{e) = i?o(I + eoffCos6l), 



(36) 



where R{9) is the radial distance from the center of the 
crystal to its interface and 9 the polar angle. The ef- 
fective anisotropy e^s represents the modification of the 
anisotropy e due to the grid. Fig. (||) illustrates a crystal 
grown to equilibrium using an input anisotropy e = 0.04. 
Using Eq. (^) we found Ccff = 0.041, within 5% of e. 
Similar accuracy was found for e — 0.02,0.03 and 0.05. 

We also tested for grid anisotropy by rotating the grid 
by 45 degrees, which should represent the lowest accu- 
racy for square elements. We compared the tip speed 
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of dendrites grown in this direction to that of dendrites 
whose principal growth direction is along the x-axis. Fig. 

shows the tip velocity for the case of a dendrite grown 
at A = 0.55 (e = 0.05, D = 2, P = 0, dt = 0.016, 
Axniin — 0.4) compared with the same case when growth 
occurs along the x-axis. The tip velocity approaches an 
asymptotic value that is within approximately 5% of the 
tip velocity computed when the anisotropy is aligned with 
the a;-direction. 



V. DENDRITIC GROWTH USING ADAPTIVE 
GRIDDING 

In this section we present results for two dimensional 
solidification with and without interface anisotropy. We 
illustrate the robustness of our algorithm and use it, in 
particular, to investigate dendritic growth at low under- 
cooling, presenting new results on dendrite tip-speed se- 
lection. 







C. Convergence and Grid Resolution 

We tested the convergence of solutions as a function of 
the minimum grid spacing Axmin- We used an undercool- 
ing of A = 0.55, with D ^2, dt = 0.016, Axmin = 0.4, 
and set A to simulate /? = 0. The parameter 7 = 1.8, 
which assured that regions of rapid change of and U 
were always encompassed in the regions of highest grid 
resolution. We examined the tip speed of a dendrite for 
0.3 < Axmin £ 1-6, finding relatively good convergence 
of the tip speed to theoretical prediction of microscopic 
solvability theory discussed above. Fig. (||) shows the 
asymptotic steady state tip velocity for each case, super- 
imposed on the solid line, which is the result of solvabil- 
ity theory for A = 0.55. It is surprising that the solu- 
tion convergence is rather good even for Axmin — 1-6. 
We have found similar convergence properties for the 
case of A = 0.25. Specifically, using Axmin = 0.4 and 
Aa;,nin = 0.78 gives essentially identical results. 

The introduction 7 in the error function ^E* gives us the 
freedom to tune the degree to which the fine grid layer- 
ing encompasses the thermal field as well as the </> field. 
Setting 7 = leads to a uniform-like mesh at the highest 
level of refinement only around the most rapidly chang- 
ing regions of (j), while the U field becomes encompassed 
in a rather disorderly combination of quadrilateral and 
triangular elements. We found that this effect can in- 
crease the tip-speed error by several percent, as well as 
increase fluctuations in tip speed. Increasing 7 produces 
a smooth layering of coarser uniform-like meshes ahead 
of the (/)-field, corresponding to region of large gradients 
in U. Fig. ^ compares the mesh around the tip of a den- 
drite grown at A = 0.65 for 7 = and 7 = 4. The figure 
illustrates the gradual mesh layering encompassing the 
thermal field for 7 = 4. In Fig. (|), D = 1, dt = 0.016, 
Axniin = 0.4 and A is chosen to simulate /3 = 0. Fig. |^ 
also shows the tip speed for A — 0.55 for the cases 7 = 
and 1.8, while Fig. (|l^) shows the tip velocity for a den- 
drite grown at A = 0.3 with 7 = and 20, respectively. 
In Fig. © £> = 10, Axmin = 0.4, dt = 0.048 and /? 0. 
In this case the higher value for 7 allows the tip velocity 
to approach within approximately 5% of the solvability 
answer, as in Ref. |32[| . Raising 7 further does not pro- 
duce any further changes in tip speed. 



A. Dendritic Growth without Surface Tension 
Anisotropy 

When the anisotropy parameter e in Eqs. (||) is set 
to zero, solidification proceeds without the emergence of 
any preferred direction. In this case it is well known that 
a seed crystal larger than a critical radius will eventually 
grow to become unstable to fluctuations, and will break 
into surface undulations via the MuUins-Sekerka insta- 
bihty 1^. Fig. |ll] shows a series of time steps in the 
evolution of a solidifying disk grown at e = 0, A = 0.65, 
D = 4, and A set to generate /3 — 0, making do = 0.1385. 
We use 11 levels of refinement and an 800 x 800 system, 
making Axmin = 0.4. For coarser meshes the Mullins- 
Sekerka instability sets in sooner due to grid noise. As 
Axmin is made smaller, grid noise becomes smaller, and 
one must wait longer for the true "thermal noise" to set 
in and make the crystal interface unstable. The dynam- 
ically evolving grids are also shown. This figure clearly 
demonstrates how the grid creation scales with the ar- 
clength of the solidifying surface. 



B. Dendritic Growth with Surface Tension 
Anisotropy 

When surface tension anisotropy is present a crystal- 
lizing disk forms dendritic branches which travel along 
the symmetry axes of the anisotropy, driven by the 
anisotropy to a steady state tip velocity []6|,p]-^ . As 
a verification of our algorithm we measured tip- velocities 
and shapes for dendrites grown at intermediate under- 
coolings. Fig. |l2| shows the dimensionless tip velocity 
(Vdo/D) versus time for A = 0.45 and 0.65 and e = 0.05. 
In Fig. (12) the dimensionless diffusivities D = 3 and 1 
and the dimensionless capillary length arc do ~ 0.186 
and 0.544, respectively. In both cases A has been set to 
simulate (3 = kinetics at the interface, while 7 = 4 
and 1.8, respectively. These values of 7 are chosen so 
as to minimize grid-layering error. The solid horizontal 
lines represent the theoretical values obtained from mi- 
croscopic solvability theory. In all cases the converged 
velocities are within a few percent of the theoretical pre- 
diction. The case of A = 0.65 includes data for three 
systems sizes. These results of system size are typical 
for intermediate A, showing a relatively rapid leveling to 
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an asymptotic speed to within a few percent of that pre- 
dicted by solvabihty theory. Fig. (|l^) also shows a plot 
of the dendrite tip shapes produced by our simulations, 
superimposed on the shapes predicted by solvability the- 
ory. 



C. Dendritic Growth at Low Undercooling 

At lower undercooling we encounter significant finite- 
size effects which cause the tip velocity to deviate from 
the solvability prediction. Fig. ^ shows the evolution 
of the tip-velocity for A = 0.25 in two different system 
sizes. For a system of size ~ 6400 x Ly = 400, the 
velocity goes to within a few percent of the solvability 
prediction. For a system size 6400 x 3200 the tip ve- 
locity seems to settle close to a value that exceeds the 
solvability prediction by 8%. This effect is even larger at 
A — 0.1, also shown in Fig. (14), where the tip speed ap- 
proaches a value about 3 times larger than that predicted 
by solvability theory. 

To understand this finite-size dependence of tip ve- 
locity at low undercooling, we note that at low A, the 
thermal fields of the two dendrite branches overlap, pro- 
ducing a thermal envelope very different from that which 
emerges for the single, isolated dendrite branch assumed 
in solvability theory. At large undercooling, each den- 
drite arm quickly outruns the other's thermal boundary 
layer, and solvability theory should apply, as is seen in 
Fig. (jj) where A = 0.7. The conditions of solvability 
theory can be better approximated at lower undercool- 
ing if simulations are performed in a domain which is 
small in one direction. For the simulation performed 
with A = 0.25 in a small box (6400 x 400), the branch in 
the y-direction is extinguished by its interaction with the 
wall and the velocity quickly approaches the solvability 
prediction. However, when both branches are present, 
as in the simulation with A = 0.25 in the larger box 
(6400 X 3200), their interaction leads to an increased tip- 
velocity because the dendrites are embedded in a circular 
rather than parabolic diffusion field. 

This is also clearly seen for a dendrite growing at 
A = 0.1 in Fig. ([T^), where the dendrite shape and its as- 
sociated field are shown for A = 0.1 {D — 13, c?o = 0.043, 
e = 0.05, Ax = 0.78, dt = 0.08). The dendrite arms 
never became free of each other in this simulation, caus- 
ing the observed deviation from solvability theory shown 
in Fig. (p^). We note that to avoid having the ther- 
mal field feel the effect of the sides of the box we per- 
form our simulations in computational domains for which 
ia; ~ (5 — 10)D/Vn- To mcct this criterion the simula- 
tion for A = 0.1 was performed in a 102400 x 51200 
domain, which is about IQD/Vn. We note that the ratio 
of the largest to smallest element size in this simulation is 
2^^. A fixed mesh having the same resolution everywhere 
would contain 9 x 10^ grid points. 

We can estimate the time t* when the growth of the 



dendrite tip crosses over from the transient regime where 
the branches interact to that where they become inde- 
pendent by equating the length of the full diffusion field, 
SiDey/"^, to the length of a dendrite arm, Vnt*. This 
gives the crossover time as 



t* = 9D/VI 



(37) 



The values for t* corresponding to the cases A = 
0.45,0.550.65, and A = 0.25 and 0.10 in Figs. (|, |l| 
and ^ confirm this scaling. 

These results at low undercooling have important im- 
plications when comparing theory to experimental obser- 
vations. In particular, since the transient time t* ^ oo 
as A ^ 0, it does not appear likely that independent pre- 
dictions for tip speed and radius, as given by solvability, 
are likely to be observed experimentally. In this regime, 
the appropriate theory to use to obtain predictions of 
the tip speed and velocity is one which explicitly takes 
into account the long range effects of interacting ther- 
mal fields of other branches. Almgren, et al. present one 
such approach |47|. In particular, study of real dendrites 
with sidebranches, growing at low undercooling, will re- 
quire such treatment. In closing, we note that while the 
independent predictions of tip speed and radius deviate 
from that of solvability theory at low undercoolings, the 
dimensionless stability parameter a* = 2doD /VnR^ does 
agree within a few percent to solvability theory. 

Further investigation of the tip speeds at low under- 
cooling, comparison with experiments and new results 
for two-sided directional solidification will be reported in 
forthcoming publications. 



VI. CONCLUSION 

In this paper we present an efficient algorithm used to 
study solidification microstructures by adaptive refine- 
ment on a finite element mesh, and solving the phase- 
field model given by Eqs. (^. Our algorithm was made 
particularly robust by using dynamic data structures and 
pointer variables to represent our evolving grid. As well, 
the modular nature of our code offers an efficient method 
of expanding the code to different situations. 

We found that our solution time scales with the ar- 
clength of the interface being simulated, allowing simu- 
lation of much larger systems and at very low undercool- 
ings. In particular, simulations for undercoolings as low 
as A = 0.1 are quite straightforward in systems larger 
than lOD/Vn- This undercooling represents the upper 
limit of dendrite growth in experiments 

Dendrite tip-velocities at intermediate to high under- 
coolings were found to agree with solvability theory to 
within a few percent. At low A, we found that the 
transient interaction of thermal fields from perpendicu- 
lar dendrite branches modifies the tip-velocity from that 
given by solvability theory at times shorter than an es- 
timated crossover time. Since this crossover time itself 
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becomes larger as A decreases, it is likely that transient 
effects will play a leading role in determining the tip ve- 
locity at low undercooling. Furthermore, this suggests 
that at low A the tip-velocity in the presence of side- 
branching will be different than that predicted by solv- 
ability theory. 

Our algorithm is currently being used to examine di- 
rectional solidification in models with unequal diffusivi- 
ties in the solid and liquid phases. These results will be 
presented in upcoming publications. 
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FIG. 1. An illustration of the quadtree element data struc- 
ture. The first frame shows an element, and four child ele- 
ments. Splitting of one of the children and one its children 
is shown, along with the branch evolution of the quadtree. 
Branches with triangles indicate square elements which are 
bridged with triangular or rectangular elements. 
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Global node numbers 
Map solution data onto work vectors 



SINGLE TIME STEP 



' ASSEMBLY AND SOLUTION 1 

Assemble finite element equations 
Solve explicitly for phase-field 
Solve implicitly for temperature 
Update solution vectors 



, i , 

-( INCREMENT TIME, M = M + 1 ) 



ERROR ESTIMATION 
Traverse quadtree and compute flux variation 
Compute error estimate on each element 



REGRIDDING 
Add/remove quadrilateral elements 
Update G-lists 

Generate new triangular and rectagular elements 



-(^ NEXT TIME SEQUENCE, N = N + T 



FIG. 3. A flow chart illustrating the algorithm program 
modules. 



FIG. 2. Illustration of all possible configurations requiring 
completion with triangular and/or rectangular elements. 




FIG. 4. A dendrite grown using the adaptive-grid method 
for A = 0.7, D = 2, e = 0.05. Clockwise, beginning at 
the upper right the figures show contours of the {/-field, the 
contour = 0, contours of the 0-field and the current mesh. 
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FIG. 5. CPU time vs. the system size, illustrating the com- 
puting time for a dendrite to move through the system of 
linear dimension Lb using our adaptive mesh method. 
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FIG. 6. The equilibrium shape of the interface, for an in- 
put anisotropy e = 0.04. The measured effective anisotropy 
e = 0.041. 
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Solvability theory 
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FIG. 9. The finite element mesh around a dendrite branch 
growing at A = 0.65, showing the grid configuration for 
(a) 7 = 4 and (b) 7 = 0. The grey-shaded lines represent 
isotherms ranging from —0.65 < U < 0.02. 



FIG. 7. The time evolution of the tip-velocity of a den- 
drite growing in the presence of surface tension anisotropy for 
A = 0.55. Data is shown for the cases where the dendrite is 
moving in the .x-dircction with two grid layering patterns, and 
along the 45 degree line. The horizontal solid line represents 
the analytic prediction of microscopic solvability. 



solvability theory 

Q □ A=0.3, ^20 
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FIG. 10. The tip-velocity of a dendrite for A = 0.3. Data 
are shown for two grid layering patterns. The horizontal solid 
lines represent the analytic prediction of microscopic solvabil- 
ity. 
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FIG. 14. The time evolution ol the tip-velocity for under- 
cooling A = 0.25 and 0.10. The data have been shifted up by 
0.025 for clarity. 



FIG. 11. The evolution of a crystal growing without intcr- 
facial anisotropy. The = contours are shown, superim- 
posed on the finite element grids. Time advances from left to 
right, top to bottom. 
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FIG. 12. The time evolution of the dimcnsionlcss tip ve- 
locity for A = 0.45 and 0.65. The horizontal lines represent 
solvability theory. The A = 0.65 case includes data for three 
system sizes. The data have been shifted up by 0.025 for 
clarity. 
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FIG. 15. Dendrite mesh and isotherms for undercooling 
AT = 0.1. (a) shows the full domain whose dimensions are 
102, 400 X 51, 200. The growing dendrite is in the lower left 
corner, (b) a close up displaying the dendrite tips, approx- 
imately 1,300 units from the origin, while the temperature 
field spreads to about 5, 000 units. 



FIG. 13. The asymptotic dendrite tip shapes for A = 0.3, 
0.45 and 0.55 (data points). The dashed lines are the shapes 
predicted by solvability theory. 
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